function [B_opt, p_opt, min_ic_val] = est_order_by_ic(x, nonzero_inds, ic_fun)
	[N,n] = size(x);

	nonzero_inds = eye(n);
	B_opt = [];
	p_opt = [];
	min_ic_val = Inf;
% TODO how to best choose p_max? just made this up
	p_max = 2*n;
	cvx_quiet(1);

	for p=0:p_max
		C = windowed_est(x, p);
		B = constrained_mlm(C, n, p, nonzero_inds);
		X = B'*B;
		ic_val = ic_fun(X, C, n, p, N)
		if ic_val < min_ic_val
			min_ic_val = ic_val;
			B_opt = B;
			p_opt = p;
		end
	end
end
